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Abstract. This article is an attempt to provide a self consistent picture, in- 
cluding existence analysis and numerical solution algorithms, of the mathematical 
problems arising from modeling photocurrent transients in Organic-polymer Solar 
Cells (OSCs). The mathematical model for OSCs consists of a system of nonlinear 
diffusion-reaction partial differential equations (PDEs) with electrostatic convec- 
tion, coupled to a kinetic ordinary differential equation (ODE). We propose a 
suitable reformulation of the model that allows us to prove the existence of a so- 
lution in both stationary and transient conditions and to better highlight the role 
of exciton dynamics in determining the device turn-on time. For the numerical 
treatment of the problem, we carry out a temporal semi-discretization using an 
implicit adaptive method, and the resulting sequence of differential subproblems 
is linearized using the Newton-Raphson method with inexact evaluation of the 
Jacobian. Then, we use exponentially fitted finite elements for the spatial dis- 
cretization, and we carry out a thorough validation of the computational model 
by extensively investigating the impact of the model parameters on photocurrent 
transient times. 

1. Introduction and Motivation 

A continuously growing attention has been paid over the last years by the in- 
ternational community and government authorities to monitoring the effect of the 
increase of global concentrations of carbon dioxide, methane and nitrous oxide on 
the quality of our everyday life. The results of the investigation carried out by the 
Intergovernmental Panel on Climate Change [1] have brought the European Union 
(EU) to the decision that carbon dioxide emissions should decrease by 20 percent, 
and that 20 percent of the energy produced in EU should originate from renewable 
energy sources, such as wind, water, biomass, and solar, not later than 2020 [2]. 
In this perspective, research and design of third generation (30) photovoltaic de- 
vices [3] for solar energy conversion into electrical and thermal energy turns out to 
be a central topic in the wider area of renewable energy sources. Roughly speak- 
ing, 3G photovoltaic devices can be divided into two main classes: electrochemical 
cells in El E] and organic-polymer cells [3 El |9] which are the topic of the present 
article. Most of investigation activity in solar cell design is devoted to the experi- 
mental study of innovative materials for efficient and flexible technologies, and is not 
presently accompanied by a systematic use of computational models to predict and 
optimize their performance. This article is an attempt to fill this gap by introducing 
the numerical engineering community to the mathematical problems that arise in 
the context of modeling and simulation of OSCs. With this aim, we try to provide a 
reasonably self-contained picture of the topic, including a discussion of the peculiar- 
ities of the model, an analysis of the existence of a solution, and the description of 
a robust computational algorithm to compute such solution. In particular, we focus 
on a special class of OSCs, namely that of Bulk Hetero-Junction (BHJ) devices, that 
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currently represent the most promising technology in terms of energy conversion ef- 
ficiency [9l E]. Charge transport in BHJs is described by a set of nonlinear PDEs of 
diffusion-reaction type with electrostatic convection coupled with a kinetic ODE for 
the temporal evolution of exciton concentration in the cell [TUl [TTl [T21 US]- Sect. |2] 
is devoted to the description of the structure and working principles of BHJs while 
in Sect. [3] the mathematical model is introduced and the connection between its 
features and the physical phenomena involved in photocurrent generation is drawn. 
Some effort is also put into highlighting the main differences between the problem 
at hand and the case of more standard crystalline inorganic semiconductor devices. 
In Sect. |4| under suitable assumptions on the model coefficients, i) we prove the 
existence of a solution of the problem in stationary conditions; and ii) we derive a 
simplified model in transient conditions, that is amenable for a qualitative analysis of 
the time response of the device, and for which we again prove existence of a solution. 
For the numerical treatment of the problem, which is the topic of Sect. |5| we carry 
out a temporal semi-discretization using an implicit adaptive method, and the re- 
sulting sequence of differential subproblems is linearized using the Newton-Raphson 
method with inexact evaluation of the Jacobian. Then, we use exponentially fitted 
finite elements for the spatial discretization, to ensure a stable approximation of the 
internal and boundary layers arising in the distribution profile of the photogener- 
ated carriers. The numerical experiments of Sect. |6]are meant, on the one hand, to 
illustrate the complex interplay among different physical phenomena determining 
the photocurrent turn-on transient time of a realistic BHJ cell in different regimes 
and, on the other hand, to characterize the range of applicability of the reduced 
model introduced in Sect. |4j In Sect. [7] we address some concluding remarks and 
indicate possible future research directions. 

2. Bulk Heterojunction Organic Solar Cells 

Before presenting the mathematical model which is the main focus of this paper, 
a schematic description of working principle of OSCs, and in particular of those with 
a BHJ structure, is in order. For more details on the subject the interested reader is 
referred to [HIIS]- The simplest possible structure for an organic-polymer based solar 
cell is depicted in Fig. [TJ two thin films composed of a conjugated organic polymer 
and of a material with high electron affinity, usually referred to as a acceptor are 
sandwiched between one transparent {e.g. indium-tin-oxide or fiuorinated tin oxide) 
and one refiecting metal contact (usually aluminum or silver). When illuminated, 
electrons in the Highest Occupied Molecular Orbital (HOMO) in the polymer are 
promoted to the Lowest Unoccupied Molecular Orbital (LUMO) thus forming an 
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Figure 2. Bulk Hetero junction OSCs. 



electron-hole pair. Such pair, which we refer to as an exciton (Fig, maj] ), in contrast 
to what is usually the case in standard inorganic semiconductors, is electrically 
neutral and has very strong binding energy (of the order of leV) with a radius in 
the sub-nanometer range. The diffusion length Xx of a moving exciton in commonly 
used polymeric materials is of the order of a few nanometers. An exciton has a 
non- negligible chance of eventually reaching the polymer /acceptor interface only 
if it was photo-generated within a distance < Ax- In case this occurs, the built- 
in chemical potential drop produced by the difference in electron affinity between 
the two materials is strong enough to stretch the exciton driving the electron and 
hole to a distance of the order of Inm thus reducing the strength of their Coulomb 
attraction. This less tightly bound electron-hole pair is referred to in the literature 
as a geminate pair (Fig. l(b)[ ) and the energy of the bond is low enough that it 
can be overcome by the electric field induced by a small voltage difference applied 
at the contacts. The newly separated electron and hole migrate, driven by electric 
field drift and diffusion forces, to the anode and cathode, respectively, where they 



are harvested thus producing a net current (Fig. 1(c)). The currently investigated 
most promising device technology to maximize the efficiency of the photogeneration 
process is the BH J cell depicted in Fig. [2] which is produced by spin-casting both 
the polymer (usually rr-P3HT or MDMO-PPV) and the acceptor (usually some 
derivative of fullerene or inorganic nanoparticles, e.g. titanium-dioxide) from a 
common solution. This process results in a highly folded structure that has the 
advantage that all photo-generated excitons eventually reach an interface, at the 
price of reducing the effective carrier mobility because of the convoluted path that 
carriers need to travel to reach the contacts. Also, from a perspective that is more 
relevant to the topic of this paper, the highly disordered structure of BHJs makes it 
difficult to characterize model parameters, as an averaging over the highly disordered 
nanostructure of the device would be required. Therefore the typical approach is to 
estimate the parameter values experimentally and resort to numerical simulations 
to properly interpret the measurement results. 



3. The Mathematical Model 

In this section we illustrate the mathematical model of the photogeneration mech- 
anisms that drive charge transport in BHJ solar cells (see [TOl |9l [TTl [121 US])- The 
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polymer /acceptor blend is represented by a homogeneous material filling a bounded 
domain Q C M"^, d > 1, with a Lipschitz boundary F = dQ divided into two disjoint 
subregions, Yd and T^, representing the interface between metal and polymer blend 
and interior artificial boundaries, respectively. We assume that meas [To) > and 
r^) n Ftv = 0, and denote by u the outward unit normal vector along F. 



(la) { dt i^^^ 

-div Jp = Gp — Rp p 



3.1. Governing Equations. Charge transport in the device is governed by the set 
of continuity equations 

On 

-div J„ = Gn- Rn n 

dp 
di 

where Qt = x (0,T), T > 0, n and p denote the electron and hole density, 
respectively. Using from now on the symbol rj to indicate either of n or p, are 
the corresponding flux densities, Grj, are the carrier generation rates, and i?^// are 
the recombination rates. As electrons are negatively charged while hole charge is 
positive, the total current density J can be expressed as J = g {Jp — Jn) where 
g > is the magnitude of the electron charge. The charge carrier flux densities are, 
in turn, each composed of an electrostatic drift term and a diffusion term 

Jn = Dn^n — fin U V(y9 

Jp = DpVp + fippVip 



(lb) < T _ n T7^^ „ inl^r, 



being the charge carrier diffusion coefficients and /i^ the carrier mobilities. The 
electrostatic potential (p satisfies the Poisson equation 

(Ic) — div{e'Vip) = q{p — n) in VLt, 

where e is the (averaged) dielectric permittivity of the blend. Notice that, as there 
are usually no dopants in organic cells, the net charge density on the right-hand-side 



in ( Ic ) is given by the carrier densities only. We denote by X the volume density of 



geminate pairs and we express its rate of change as 

dX 

(Id) — = g-r inVtj 



The geminate-pair generation rate g in [Id] can be split into two contributions as 

(2) g = G{x,t)+Tpn, 

(a) accounting for the rate at which excitons reach the material interfaces and are 
partially separated and [b) accounting for the rate at which free electrons and holes 
are attracted to each other and recombine. Process (6) is referred to as bimolecular 
recombination and the coefficient 7 is described according to the Langevin theory [7] . 
The rate of process (a) is equal to the rate G{x,t) at which photons are absorbed, 
which we assume in what follows to be a known function of position and time. As in 
BHJ all excitons are eventually transformed in geminate pairs it is legitimate, with 
a slight abuse of notation, to use in the following the two terms as synonyms. As 



for the term r in (Id) it can also be split into two contributions as 

(3) r = kdiss-^ "I" kj-ec-^ y 
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(c) accounting for the rate at which geminate pairs that are not spht recombine 
and (d) accounting for the rate at which free electrons and holes are produced by 
separation of a bound pair. We assume the coefficient krec to be a given constant 
while kdiss depends on the magnitude of the electric field E = — Vy? as described 
in [7j. As we assume free carriers to be generated only by dissociation of a geminate 
pair and to be annihilated only by recombination into a geminate pair, the generation 
rates satisfy Gn = Gp = k^iasX while for the recombination rates = RpP = •ypn 
holds. 

We wish at this point to stress some peculiarities of the model we have introduced 
compared to the standard case of crystalline inorganic semiconductors. The main 
difference is represented by the strong infiuence that the exciton reaction kinetics de- 
scribed by equation (Id) has on device performance. Indeed, such a kinetics affects 
both the energy conversion efficiency in the steady state operation and the turn-on 
transient time. This latter, in particular, is relevant for the characterization of ma- 
terial properties that can not be determined by first-principles because of the highly 
convoluted device nanostructure. Furthermore, although equations (la)-(lb) are 
analogous to those describing charge transport in ordered inorganic semiconductors, 
the physical driving mechanisms at the microscopic level are quite different. In par- 
ticular, while in monocrystalline semiconductors charge carriers are essentially free 
to move within delocalized orbitals, in the materials we study here transport hap- 
pens via hopping of charges between localized orbitals. This microscopic difference 
is reflected in the macroscopic models for the diffusion and mobility coefficients for 
organic semiconductor materials which (i) introduce very different dependencies on 
temperature and electric fleld magnitude [HI [15], and (ii) introduce a dependency 
of the mobility on the carrier densities 



3.2. Boundary and Initial Conditions. A delicate and important issue is that of 
devising a set of boundary conditions to accurately describe the complex phenomena 
of charge injection and recombination occurring at the interface separating the 
metal contacts from the semiconductor bulk. Precisely, according to [ISl [ID] , such 
conditions are expressed in the following Robin-type form 

(4a) KnJn-l' = Pn-ann OB. T D X (0, T) 

(4b) KpJp-i/ = Pp-apP on X (0, T), 

where are non negative parameters while are the rates at which charges are 
injected into the device and a^r] are the rates at which electrons and holes recombine 
with their image charges at the contacts, respectively. Reliable models for the 
above parameters are still subject of extensive investigation as the basic description 
proposed in the milestone reference [16] needs to be modifled via empirical fltting to 
avoid the occurrence of unphysical behavior in the computed solution [121 IT]- As 
for the electric potential, the Dirichlet condition 

(4c) <^ = ^fl onr^x(0,T) 

is enforced, where the datum accounts for both the externally applied voltage 
and the work-function difference between the contact materials. On T^, which 
represents the interior artiflcial boundary, homogeneous Neumann conditions for the 
flux densities and the electric fleld are imposed. Finally, positive initial conditions 
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n{x,0) = riQ^x), p{x,0) = Po{x), and X{x,0) = Xq{x) are needed to complete the 
mathematical model. 



4. System Analysis of the Model 
In this section, we deal with the analysis of the existence of a solution of sys- 



tem (la)-(ld) in both stationary and transient regimes, under the following as- 
sumptions: 

(HI): 7, kdiss, krec and G are all positive constant quantities in Qt] 

(H2): = Vthfirj, Vth being the thermal voltage and /x^ > /U^q > a.e. in Qt] 

(H3): Vn,Vp < f"^"^ < +00 where := Hn\E\; 

(H4): = and a^, are functions of position only in (4a)- (4b). 

Although the purpose of the set of hypotheses (H1)-(H4) is mainly to reduce 
the mathematical complexity of the problem, we wish here to comment about their 
physical plausibility. Assumption (HI) allows us to express in an easy manner the 
dependent variable X as a function of n, p and of the input data G and Xq, in such 
a way that the resulting equivalent system (in the reduced set of unknowns (f, n 
and p) can be written in the form of a two-carrier drift-diffusion (DD) model. As 
the coefficients involved in (HI) depend, in general, only on the magnitude of the 
electric field, such an assumption is reasonable if the field itself varies weakly within 
the simulation domain, which is often the case in realistic photovoltaic devices as is 
confirmed by the numerical experiments of Sect, [g] Assumption (H2) is the classical 
Einstein relation valid in inorganic semiconductors and corresponds to neglecting 
the (higher order) effect of energetic disorder ^3]. The saturation of convective 
velocities expressed by assumption (H3) is reasonable in a structure that is highly 
folded as that of BHJs and is indeed commonly employed in commercial packages for 
organic semiconductor simulation [18]. Assumption (H4) corresponds to an infinite 
carrier recombination rate at the contacts. 



4.1. Stationary Regime. Setting dX/dt = in (Id), we can eliminate the depen 



dent variable X in favor of n, p and of the input function G, to obtain 

(5) X{x) = tG + '^Tp{x)n{x) 
where 

(6) r ^ 



l^diss ~l~ ^rec 

is the time of response of the generation/recombination terms to light stimuli. Us- 
ing ([5])-(|6]) and (H4), the stationary OSC model reads: 

{— div(£:V(y9) = qijp — n) 
-divJ„ = r {kdiss G --f Kec pn) in VL 

-div Jp = r {kdiss G - 7 Kec pn) , 

supplied with the boundary conditions 



(8) 



Lp = '^D, n = nD-=—, P = Pd-=— on r 



D 



p 



Jn ■ 1^ = Jp ■ 1^ = ■1^ = on r^r. 
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Theorem 1 (Existence of a solution in stationary regime). Let assumptions (Hl)- 
(H4) be satisfied and D^riD^pD) G {L°^{Vd)Y. Then, problem ((7])-(|8| admits a 
weak solution u* , v*) G {H^iVL) nL°°(r2))^ and there exist positive constants M. , 
A4, JC, )C such that 



(9) 



M <n*, p* <M, ic <ip* <k: 



a.e. in Q. 



The proof of Theorem [T] foUows closely the guidelines of Sect. 3. 3 and is 
sketched below. Using (H2) we can write the two flux densities as 



(10) 



Jp = fipVthnre-'^/^^'^Vv, 



where the new (dimensionless) dependent variables u and v are related to the carrier 
densities n and p by the Maxwell-Boltzmann statistics 



(11) n = Ur ue'*'' ""■ , p = Ur ve 

Hr > being a reference concentration. System ([7])-([8]) then becomes: 



(12) 

and 
(13) 



-divifinVth e^/^*'^Vn) 



q Uriue'^/^^'^ - t;e-'^/^"') 
[1 — uv) 



Ur 



in VL 



(1 — uv) 



^J, Pd ^n/Vfh 
(/3 = U = Ud '■= e o/yth^ v=yj^-= e 'D/''th 

Jn ■ 1^ = Jp ■ 1^ = ■ 1/ = 



on Yd 
on Ftv. 



Using the boundedness of the Dirichlet data, the positivity of ra^ and po and choos- 
ing Ur in such a way that {■ykrec'nf.) / (kdissG) = 1, we can see that 



(14) 
where 



a.e. on To, 



and 



vE^"^ := max < max (sup (-</?„!)), sup(v?pz))), - min(inf (-</?„£,), inf (v9p£))) 



(PnD ■= - Vthln^no/nr), (PpD ■= *D + Vthln{pD/nr). 



Then, by applying Theorem 3.3.16 of [19] to system (12)-([13j) and going back to the 
original variables n and p via the inversion of (11), we conclude that Theorem [l] 
holds with 



(15a) 
(15b) 

where 



K, = Ur e 



th 



th 



M. = min I inf "^d, 



M = max ( sup\E'£i, 



sup I^^dI 
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4.2. Transient Regime. Analogously to what we have done in Sect. 4.1 in the 
stationary case, we can use (Id) to eliminate the dependent variable X in favor of 
n, p and of the input functions G and Xq, to obtain 

(16) X{x,t) = ^{x,t)+-f [ p{x,s)n{x,s)e-'-^-''>^^ ds, 

Jo 

where ^{x,t) := Xo(a;)e~^/^ + r G{1 — e~^/^). The quadratic convolution term 
in ( 16 ) makes the dependence of the current on the electron and hole densities non- 
local in time with a "memory window" of size proportional to r. For the subsequent 
existence analysis it is convenient to manipulate such term so that we can write the 
continuity equations in the following equivalent form: 

- div J„ = kdiss^ -IT [krec + kdiss^i'^''') pn+ I 



(17) 

where 
(18) 




div Jp = kdissC - 7 ^ {krec + kdisse *''^) pn+ I 



I := 'yk, 



diss 



[p{x, s)n{x, s) — p{x,t)n{x,t)]e ^^^'^ ds. 



Although I is no more a convolution integral, it has the interesting property of 
vanishing both at t = and t = +oo, from which we expect, at least formally, 
that replacing the integral J by a suitable approximation, say /, should not have a 
significant impact on the model behaviour as long as it preserves the asymptotics of 
I. Our choice is to use a trapezoidal quadrature rule, yielding 

t 



(19) 



/ ^ / = 7 kdiss T^e [p{x, 0)n{x, 0) - p{x, t)n{x, t)] . 



It is easy to see that / vanishes both at t = and t = +oo; moreover, the approxi- 
mate formula ( 19 ) as the advantage of lumping the non-locality of / into a quadratic 
term that has the same form as the generation/recombination rates already present 
in the right-hand side of (17) The resulting reduced model reads: 

q{p — n) 
div Jn = Gn - RnJl 



(20) 



On 

di 
dp 



L dt - ^^"-^^ 



in Qt 

Gp — RpP, 

where the modified generation/recombination mechanisms are defined as 

Gn 



(21) 



t 



Gp = kdissi{x, t) + 7 kdiss^-e '^^^p{x, 0)n{x, 0) 



RnU = RpP = 7 



'^ij^rec ~l~ ^diss^ ^ ^ ~^ ^diss 



-t/r 



p{x, t)n{x, t). 



Having derived a new, simplified model, it is natural to ask to which extent the novel 
formulation is capable to describe correctly the main features of the performance 
of an OSC. With this aim, we first investigate the quality of the approximation 
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provided by /; the quadrature error associated with the use of the trapezoidal rule 



in ( 19 ) is given by the following relation [20j 



f22) 



12 



(*-^)/MA"(C) + -A'(C) + ^(A(C)-A(t)) 



r 



SI n( 



where ( G (0,t) and A(s) 
ligible as t — )■ or t — 7- +oo, as expected, meaning that the predicted (computed) 



s). Eq. (22) shows that E{t) becomes neg- 



stationary current is independent of the use of (19) or the exact expression (17), as 



numerically verified in Sect. |6.2[ However, for a finite value of time t, the discrepancy 
between the exact convolution term and its approximation may be non-negligible. 
A reasonable estimate of the error would require a knowledge on the temporal be- 
havior of the photogenerated carrier densities n and p as a function of time. This 
knowledge not being available, we can still gain some information on the quadra- 
ture error by an analogy with the approximation of the recombination/generation 
term that is usually carried out in the study of currents in a p — n junction in 
the inorganic case (see [2T]). This analogy suggests that the value of E(t) during 
the photocurrent transient (i.e., for t sufficiently far from but also sufficiently far 
from stationary conditions) might become significant if the OSC is operating under 
high injection conditions, or, equivalently, high current level conditions. Again, this 



latter statement is numerically verified in Sect. 6.2 



Theorem 2 (Existence of a solution in the transient regime). Let assumptions 
(H1)-(H4) be satisfied, and the initial data U := (no,Po); -^o o,'>T'd the function \1/ be 
such that U e {H^inr) n L'^^Qt))^, with U > 0, Xq e L°°{n) with Xq > 0, and 
^ G H^i^lr) nL°°(nr). Then, setting u := {n,p), system supplied with 



the initial/boundary conditions (4a)-(4c), admits a weak solution {(p,u) such that: 
(1) u > a.e. 



in 



■T, 

(2) u{x, 0) = U{x, 0) andu-U e (0, T; H^f ; 

(3) w G (C(0,T;L2(fi)) nL°°(fiT))'; 

(4) ^GL2(0,T;/7;)^- 

(5) ^ - ^ G L2 (0, T; Hq) with G (fi^) , 

where Hq := {v & H^{Q) : v\r^j = 0} and H'q is its dual. 

Moreover, using (16) and the regularity of Xq, n andp, we have that 

dX 

x, — eCio,T;L'in))nL^inT) 

with X{x,t) > for all t > and for a.e. x E Q. 

The proof of Theorem |2] consists of verifying that all of the assumptions (Ei) 
(Eiv) of [22], p. 



296 are satisfied. It is immediate to see that the functions Rr/ 
are positive for p > and n > and satisfy locally Lipschitz conditions, with a 
Lipschitz constant which is uniform in time and equal to 27. As a matter of fact, 
for any n',p', n",p" and for any x and t, we have 

\Rn{x,t,n',p') - Rn{x,t,n",p")\ < 7 {r^krec + kdiss) + rkdiss) \p' - p"\ 

< 27|p'-p"|, 



and the same estimate holds for Rp provided to exchange \p' — p"\ with 



\n 



n 



Moreover, (H2) and (H3) ensure that (20)2,3 are uniformly elliptic with uniformly 
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bounded convective velocities. Then, by applying Theorem 2 of ^22J, we conclude 
that Theorem [2] (of the present article!) holds. 



5. Numerical Discretization 

In this section, we illustrate the numerical techniques for the simulation of the full 
model ([T])-(|4]), as the same approach can be used, with slight modifications, to treat 



the reduced approximate model (20)-(21). In designing the algorithm presented 



here, our aim is twofold: on the one hand, it seems natural to try to adapt methods 
that are known to work efficiently and reliably for transient simulation of inorganic 
semiconductor devices (see, e.g., [23J Chapt. 6, Sect. 4); on the other hand, as the em- 
phasis of the present paper is on accurately estimating photocurrent transient times, 
it is necessary to apply advanced time-step control techniques |211[2S]- To this end, 
our chosen approach is based on Rothe's method (also known as method of horizon- 
tal lines) which consists of three main steps: first, the time dependent problem is 
transformed into a sequence of stationary differential problems by approximating the 
time derivatives by a suitable difference formula; then, the resulting nonlinear prob- 
lems are linearized by an appropriate functional iteration scheme; and, finally, the 
linear differential problems obtained are solved numerically using a Galerkin-Finite 



Element Method (G-FEM) for the spatial discretization. Sects. 5.1, 5.2 and 5.3 



below discuss in more detail each of these steps; it is worth noting that, with minor 



modifications, the linearization techniques of Sect. 5.2 can also be applied to treat 
the stationary model ([T]). 



5.1. Time Discretization. To transform the time dependent problem ([T])-(|4]) into 
a sequence of stationary problems, we replace the partial time derivative with a 
suitable finite difference approximation, specifically, the Backward Differencing For- 
mulas (BDF) of order m < 5 (see, e.g., [ll]. Sect. 10.1.2). To describe the resulting 
stationary problem, let = to < • • • < tx-i < tx < T he a strictly increasing, not 
necessarily uniformly spaced, finite sequence of time levels and assume the quantities 
Ml = n, U2 = p, X and ip to be known functions of x for every tj^, k = . . . K — 1. 
Then we obtain: 



(23) 



-div{eVipK) + q {uk - Pk) 

m 

^dkUK-k - divJ„(nK; VifK) - Uk 

m 

y^dkPK-k - dwJp{pK; Vcpk) - Uk 

k=0 

m 

''^OkXx-k — Wk 








0, 



k=0 
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where fk = f{x,tk) for any generic function / = f{x,t), and 
Uk := U(yipK,nK,PK,XK,tK) 

= GpiVifK, nx, Pk, Xk, tx) - RpiVifx, ^x, Pk, Xk, tx) Pk, 



K 



-- W{VipK, UK, Pk, Xk, Ik) 

g(yipK,nK,PK,XK,tK) - r(yipK,nK,PK,XK,tK). 



System (23), together with the constitutive relations for the fluxes given in (lb) 
and the set of boundary conditions (|4]), constitutes a system of nonlinear elliptic 
differential equations (la) coupled to an algebraic constraint equation (Ic). In our 



implementation, the selection of the next time level tx and of the formula's order 
m, as well as the computation of the corresponding coefficients 6k, k = 0, ...,m, 
is performed adaptively to minimize the time discretization error while minimizing 
the total number of time steps via the DAE solver software library DASPK [26l [27] . 
Notice that, if m = 1, we have 9q = —9i = ^ , ^/c = 0, k > 1, and the temporal 

semi-discretization of system ([T|-(|4]) coincides with the Backward Euler method. 

5.2. Linearization. To ease the notation, throughout this section the subscripts 
denoting the current time level will be dropped. Let y := [ip, n, p, X]^ denote 
the vector of dependent variables and let denote the null vector in M^. Then, the 



nonlinear system ( 23 ) can be written in compact form as 



(24) 



F{y) = 0, with 



fv{^,n,p) 
fn{v,n,p,X) 
fp{ip,n,p,X) 
fx{'f,n,p,X) 



The adopted functional iteration technique for the linearization and successive solu- 



tion of problem (23) is the Newton- Raphson method. One step of this scheme can 

be written as 

(25) 





dniU) 


dpiU) 











-fv{^,n,p) 




dnifn) 


dpUn) 


dxifn) 




An 




-fn{^,n,p,X) 


dMp) 


dnifp) 


dpifp) 


Mfp) 




Ap 




-fp{^,n,p,X) 


dMx) 


dnifx) 


dpifx) 


dxifx) 


{ip,n,p,X) 


AX 




_ -fx{^,n,p,X) 



where da{f) denotes the Frechet derivative of the nonlinear operator / with respect 



to the function a. More concisely, we can express (25) in matrix form as 



J{y) Ay = -F{y), 

where J is the Jacobian matrix and Ay := [Aip, An, Ap, AX]'^ is the unknown 
increment vector. The exact computation of all the derivatives in the Jacobian on 



the left hand side in (25) can become quite complicated if the full model for all 
the coefficients (most notably the electric field dependence of k^iss, fJ'n and /ip) is 
taken into account. Moreover, this would require cumbersome modifications to the 
solver code whenever a new coefficient model is to be implemented. One alternative 
could be to employ a staggered solution algorithm, often referred to as Gummel- 
type approach in the semiconductor simulation context [281 129] • The decoupled 
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approach is well known to be more robust as compared to the fully coupled Newton 
approach (25) with respect to the choice the initial guess and also less memory 
consuming. As in this particular study we can rely on the knowledge of the system 
variables at previous time levels to construct a reasonable initial guess and as we 
are dealing with an intrinsically one-dimensional problem (see Sect. [6]), memory 
occupation is not likely to be a stringent constraint, so that we adopt a quasi- Newton 
method where, rather than the exact Jacobian J{y), we use an approximation J{y) 
in which the dependence of the mobilities, of the diffusion coefficients and of the 
dissociation coefficient on the electric field is neglected. This approach has the 
further advantage of facilitating the use of a standard software library like DASPK 
for advancing in time. 



5.3. Spatial Discretization and Balancing of the Linear System. Once the 
linearization described in the previous section is applied, the resulting linear sys- 
tem of PDEs is numerically approximated by means of a suitable G-FEM. Pre- 
cisely, to avoid instabilities and spurious oscillations that may arise when the drift 
terms become dominant, we employ an exponential fitting finite element discretiza- 
tion [30l [3ll [32l |33] . This formulation provides a natural multidimensional extension 
of the classical Scharfetter-Gummel difference scheme [3H [35] and ensures, when ap- 
plied to a carrier continuity equation in the DD model, that the computed carrier 
concentration is strictly positive under the condition that the triangulation of the 
domain Q is of Delaunay type. It is important to notice that, when implementing 
on the computer the above described procedure, the different physical nature of the 
unknowns of the system and their wide range of variation may lead to badly scaled 
and therefore ill-conditioned linear algebraic problems, which in turn can negatively 
affect the accuracy and efficiency of the algorithm. To work around this issue, we in- 



troduce two sets of scaling coefficients, denoted {a^, 
and restate problem (24) as 



o-n, o-p, cTx} and {9?, n, p, X}, 



(26) 



(Te- 



ar, 



-fni'^0,nn,pp,XX) = 
-fp{(p(p,nn,pp,XX) = 



fx{(p0,nn,pp,XX) = 0, 



where (p := (p/(p, n := ra/n, p := p/p and X := X/X. Solving (26) for the scaled 
dependent variables [if, h, p, X]'^ corresponds to solving a system equivalent to (25) 



where the rows of the Jacobian J and of the residual F are multiplied by the factors 
{l/cr<^, l/cr„, l/cTp, l/cTx} while the columns of J are multiplied by the factors 
{(^, n, p, X}. Computational experience reveals that a proper choice of the scaling 
coefficients might have a strong impact on the performance of the algorithm. For 
example, to obtain the results of Fig. |4] a suitable choice was found to be that of 



settmg 



(Jr. 



10^ (Tx = 10^ and ^ = 1, 



n 



10^^ X 



10 
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while values differing by more than one order of magnitude from such choice were 
found to hinder the ability of the DAE solver to reach convergence. 
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6. Numerical Results 



This section is devoted to presenting the results of numerical simulations carried 
out with the algorithms described in Sect. |5| In particular, in Sect. 6A_ we discuss 



the simulation results for a realistic BHJ device focusing on the impact of the model 
parameter values on the turn-on transient time in different operation conditions. 



Sect. 6.2 is devoted to characterizing the region in the model parameter space where 



the approximate formula (19) and the resulting reduced model (20) are reliable. In 
both cases the considered device has a thickness Lqsc = 70nm and the contact 
materials are ITO and Al for the transparent and reflecting contact, respectively. 
As no external voltage is applied to the device, this results in a total voltage drop 
across the device AV = 0.5V. The relative permittivity constant is = 4 and the 
operating temperature is 3007^. As the thickness of the device is much smaller with 
respect to the dimensions in the other directions (typically many orders of magnitude 
larger) and the donor /acceptor blend is considered to be uniform, the simulations 
presented here are performed in one spatial dimension, so that the computational 
domain is modeled as the segment Q = [0,Losc] with the cathode at a; = and 
the anode at x = Lqsc- Also, as the device length is quite small compared to the 
wavelength of visible light, it is reasonable to consider the photon absorption rate 
G to be constant in Q at any t G [0, T]. 



6.1. Simulation of a realistic device. In this section we present simulation re- 
sults of the realistic BHJ device whose data are given in |12] . The computations try 
to reproduce the measurements that are commonly performed in research laborato- 
ries to characterize the device material properties and are meant to show the ability 
of the model to capture the complex dependence of the turn-on transient time on 
both the mobility coefficients and the exciton dissociation/recombination dynamics, 
and the predominance of one or the other of such phenomena depending on the op- 
eration conditions, i.e. on the intensity of the light to which the device is exposed. 
Throughout this section we use for the coefficients in the boundary conditions (|4]) 
the current injection model of [161 110] corrected as in fl2[ [T7] to increase the carrier 
surface recombination rate, thus avoiding the occurrence of spurious charge build- 
up effects near the contacts. The exciton dissociation coefficient kdiss is considered 
to depend on the electric field according to an Onsager-like model given by the 
nonlinear formula presented in [lOj with the initial separation of the geminate pair 
set to a = 1.5nm, while the recombination krec rate is constant. The bimolecular 
recombination coefficient 7 depends on the carrier mobilities and on the material 
permittivity e as resulting from Langevin theory \7], therefore, as we consider here 
the carrier mobilities to be constant, 7 is a constant as well. Figure [3] shows the 
photocurrent evolution in response to an abrupt turn-on of a light source; for each 
row in the figure the charge carrier mobilities are kept constant while the exciton 
recombination coefficient is varied whereas for each row in the figure the mobilities 



are fixed and the recombination coefficients vary. By comparing Figs. 3(a) and 3(c 



to Figs. 3(b) and 3(d) one can notice that the strong impact of the recombination 
rate coefficient krec on the transient duration in high illumination conditions (dashed 
lines) completely overshadows the effect of the carrier transport properties, while in 
low illumination conditions (solid lines) the importance of the effect of krec is less 
apparent so that the transient time is more related to the value of the mobilities. 
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Figure 3. Transient currents at low and high intensities with differ- 
ent mobihties and exciton recombination rate coefficients. For (a) and 
(b) the mobihty was 2 x 10~^cm^V^~^s^^ with geminate recombination 
rate constants k^ec = 1 x 10^s~^ and 1 x 10^s~^ respectively. For (c) 
and (d) the mobility is 2 x IQ'^crn^V^^ s^^ with krec = 1 x lO'^s"^ and 
krec = 1 X 10^s~^ respectively. 



Figure |4] shows the time evolution of the electron density in the device under 
strong illumination conditions {G = 4.3 ■ 10^° m~^s~^). Hole density is not shown 
in the figures because, due to the choice of equal mobilities, it is the exact mirror 
image of the electron density. As previously mentioned, due to the absence of fixed 
charges (dopants) within the bulk of the device the charge densities do not show the 
steep interior layers that are the main peculiarity of inorganic semiconductor models 
and lead to the main difficulties in the numerical simulation of such devices. Also 
the steepness of the boundary layers is less extreme in the case of organic devices 
and is further mitigated by the inclusion of finite surface recombination speed in the 
boundary conditions. 

The consistency of the results shown here with those of [12j is a strong indication 
of the robustness of the numerical algorithm of Sect. [5j Finally, Figure [5] shows 
the magnitude of the electric field along the device, for low illumination (solid line) 
the electric field is practically constant throughout the device while for high light 
intensity (dashed line) its deviation around its mean value {E) = AV/Lqsc is about 
30%. 
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Figure 4. Time evolution of the electron distribution at high in- 
tensity with (a) high charge generation efficiency and (b) low charge 
generation efficiency. 
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Figure 5. Value of the computed electric field for a device with mo- 
bilities fiy = 2 X 10^^cm'^V~^s~^ and recombination rate constants 
k-ppp — X X X tS 



6.2. Validation of the simplified model. In this section we wish to estimate 
the impact of the approximation (19) on the simulation results for parameter values 
within a physically plausible range. To be consistent with assumptions (H1)-(H4) 
of Sect. [4], throughout the present section we enforce that all model coefficients be 
constant by replacing the spatially varying electric field E in the coefficient models by 
its mean value (E) = —AV/Losc- Furthermore we consider carrier recombination 
at the contacts be instantaneous, so that the boundary conditions Q degenerate 
into simple Dirichlet type conditions. The plausibility of these assumptions has been 
already addressed at the beginning of Sect. |4] and in the discussion of the numerical 
results of Sect. O In all subsequent figures, the dashed line refers to the solution 
computed with the full (3 carrier) model dl])-rt4| while the solid line refers to the 
simplified approximate (2 carrier) model (20)-(21). 

Figures [6]j7jj8] refer to a device under low light intensity conditions and show the 
impact on the turn-on transient time of the value of the mobilities, of the geminate 
pair dissociation rate and of the recombination rate, respectively. 
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Figure 6. Photocurrent transient at low light intensity: effect of 
mobility on rise time. 




Figure 7. Photocurrent transient at low light intensity: effect of 
dissociation rate on rise time. 
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(a) krec = lO^S-l (b) krec = lO^fi-l 



Figure 8. Photocurrent transient at low light intensity: effect of 
geminate pair recombination rate on rise time. 

One may observe that, while at low intensity a change of one order of magnitude 
in the value of the mobility produces an almost equal change in the transient time, at 
high light intensity (Fig. [9]) a similar change in the mobility has an almost negligible 
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impact. In this latter regime, variations in the dissociation rate kdiss (Fig- 10) and, 
more notably the recombination rate krec (Fig- H), produce a more dramatic effect. 
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Figure 9. Photocurrent transient at high light intensity: effect of 
mobility on rise time. 
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Figure 10. Photocurrent transient at high light intensity: effect of 
dissociation rate on rise time. 




Figure 11. Photocurrent transient at high light intensity: effect of 
geminate pair recombination rate on rise time. 
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The analysis of the above resuhs displays the complex relation between the tran- 
sient behaviour of the device and the strongly nonlinear interplay among the several 
occurring physical phenomena and shows the ability of the simplified model (20)- 
(21 ) to capture such behaviour in most circumstances. The only situation where the 
two models disagree is in the case of a device with high generation efficiency (i.e., a 
low value of krec) under high light intensity (cf. Fig. [lT|a)). Finally the steady-state 
current predicted by the reduced model is always in perfect agreement with that of 
the full model, as expected. 



7. Conclusions and Future Work 

In this article, we have dealt with the mathematical modeling and numerical sim- 
ulation of photocurrent transients in nanostructured mono-layer OSCs. The model 
consists of a system of nonlinear diffusion-reaction PDEs with electrostatic convec- 
tion, coupled to a kinetic ODE. We have proposed a suitable reformulation of the 
model which makes it similar to the drift-diffusion system for inorganic semiconduc- 
tor devices. This has allowed us to prove the existence of a solution for the problem 
in both stationary and transient conditions and to highlight the role of exciton dy- 
namics in determining the device turn-on time. For the numerical treatment, we 
carried out a temporal semi-discretization using an implicit adaptive method, and 
the resulting sequence of differential subproblems was linearized using the Newton- 
Raphson method with inexact Jacobian. Exponentially fitted finite elements were 
used for spatial discretization, and a thorough validation of the computational model 
was carried out by extensively investigating the impact of the model parameters on 
photocurrent transient times. 

Future work is warranted in the following three main areas: 1) extensions to 
the model; 2) improvement of the analytical results; and 3) development of more 
specialized numerical algorithms. In detail: 

1) : we intend to include exciton transport in order to be able to simulate 
multi-layer or nanostructured devices [HI |Hl EHl ET] ; 

2) : we aim to extend Theorem |2] to cover the full problem ([T|-(|4]). A possible 
approach to achieve this result is to apply Theorem |2] locally on a partition 
of [0, T] into sub-intervals of size At, and verify the hypotheses of the Aubin 
lemma [38] to extract a limiting solution as At —t- 0; 

3) : starting from the above idea, we intend to devise a numerical algorithm for 
the local approximation of the full model system over each sub-interval of size 



At using the reduced model (20)-(21 ). The computer implementation of this 
approach is straightforward as it basically amounts to a successive applica- 
tion of the formulation discussed in Sect. [5] on each time slab. Furthermore, 
we intend to improve the robustness of the nonlinear solver with respect 



to the choice of scaling parameters (cf. Sect. 5.3) by adopting a staggered 
solution scheme based on some variant of Gummel's Map [281 ESj. Such 
scheme could be either employed as an alternative to the current Newton 
solver or, even more effectively, combined with this latter in a predictor- 
corrector fashion. The above modifications to the solution algorithm are 
of great importance in dealing with the simulation of the multidimensional 
device structures mentioned at item 1). 
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